function handles = excitationProfileAnalysis(handles, spacingXY, somaXYZ, figFlag)
% excitationProfileAnalysis     Analyze responses in terms of distance from soma
% Syntax:
%   excitationProfileAnalysis(maps, spacingXY, somaXYZ, figFlag)
% - maps: a map stack.
% - spacingXY is a 2-element vecor
% - somaXYZ is a 3-element vectors.
% - NB -- somaXYZ(3) is the focal plane relative to the soma plane
% - set figFlag to 1 for graphical display of map and gaussian fit
%
% Directly derived from responseDistance.m, with addition
% of stuff from analyzeArray
%
% Private version for mapAnalysis3p0
%
% See also responseDistance, analyzeArray

% note -- need to figure out how to analyze correctly for Z > 0
% correct way is probably
% - do the gaussian fit without factoring in Z
% - factor the Z in for the mean-distance calculation

% more notes -- 
% - could re-write distance arrays using meshgrid
% - need to revise distance arrays to accommodate non-zero soma positions

% gs Sept 2004 -- Private version for mapAnalysis3p0
% gs may 2006 -- modified for compatibility with ephus software
% GS20060530 - fixed error involving distance calculations (relating to EP stacks for 2P LSPS)
% ********************************************************************
% global X Y

disp(' ');
disp('% EXCITATION PROFILE ANALYSIS:');
disp(' ');

% if nargin == 0 
%     data = getData('EP_adult_rat_L4_092404'); 
%     maps = data.maps;
%     spacingXY = data.spacingXY;
%     somaXYZ = data.somaXYZ;
%     figFlag = 1;
% end

maps = handles.data.map.mapActive.mapDataArray(:,:,6);
[R,C,P] = size(maps);

% DISTANCES ----------------------------------------------------------

% arrays and vectors of distances
colVals = [-((R-1)/2)*spacingXY(1):spacingXY(1):((R-1)/2)*spacingXY(1)]';
colVals = colVals - somaXYZ(1);
colArray = repmat(colVals, 1, C);

rowVals = [-((C-1)/2)*spacingXY(2):spacingXY(2):((C-1)/2)*spacingXY(2)];
rowVals = rowVals - somaXYZ(2);
rowArray = repmat(rowVals, R, 1);

% planeArray = somaXYZ(3) * ones(R,C); % GS20060531 - commented; no longer needed (& causing error)

% distArray = sqrt(rowArray.^2 + colArray.^2 + planeArray.^2); % GS20060530 - commented; last term was wrong
distArray = sqrt(rowArray.^2 + colArray.^2);
% dists = round(distArray), return
[r,c] = size(distArray);
distances = reshape(distArray, r*c, []);
distances2 = [-distances; distances];

% RESPONSES ----------------------------------------------------------

for p = 1:P
    map = maps(:,:,p);
    [r,c] = size(map);
    responses = reshape(map, r*c, []);

    % data format for later plotting and gaussian fit
    responses1(:,p) = responses;
    responses2(:,p) = [responses; responses];

    % distances of responses > 0
    I = find(responses > 0); 
    
    % display raw data
    disp('#   R');
    disp('------');
    disp(num2str(round([responses(I) distances(I)])));
    disp('------');
    
    % mean distance weighted by response number
    resps = responses(I);
    dists = distances(I);
    distVec = [];
    for k = 1:length(I)
        distVec = [distVec; repmat(dists(k), resps(k), 1)];
    end
    
    meanWeightedDistance(p) = mean(distVec);
    stdWeightedDistance(p) = std(distVec);
    numel(distVec);
    semWeightedDistance(p) = std(distVec) / sqrt(numel(distVec));
    
    % unweighted version
    meanDistance(p) = mean(distances(I));
    stdDistance(p) = std(distances(I));
    semDistance(p) = std(distances(I)) / sqrt(numel(I));
end


% response distance analysis in terms of each unique distance

% need to re-do this -- if one map, error bars are within-map SEM; if >1 map, bars are across-cell SEM

U = unique(distances)'; % list of unique distance values
for n = 1:numel(U)     % for each unique distance ...
    I = find(distances==U(n));   % get the indices for occurrences of that distance
    % distances(I) % uncomment to verify correct indices
    A = [];
    A = responses1(I, :);
    A = reshape(A, [], 1);
%     Aavg(n) = sum(A) / P; 
    Aavg(n) = mean(A);
    Asem(n) = std(A)/sqrt(numel(A));
end
Xuniq = U;



% average map
mapAvg = mapAverage(maps, 1);

if size(maps, 3) == 1 % for now, only if 1 EP sent to function
    p = 1;
    % analysis
    totalNumberOfSites = length(find(maps(:,:,p)));
    totalNumberOfSpikes = sum(sum(maps(:,:,p)));
    normNumber = totalNumberOfSpikes * (spacingXY(1)/1000) * (spacingXY(2)/1000);
    %     spikesPerSite = round(100*sum(state.uncaging.map.apNumber(state.uncaging.map.apNumber>0)')/length(tracesWithAPs'))/100;
    %     meanDist = round(mean(state.uncaging.map.distFromSoma(state.uncaging.map.apBinary')));
    warning off MATLAB:divideByZero;
    spikesPerSite = totalNumberOfSpikes / totalNumberOfSites;
    meanDist = round(meanWeightedDistance(p));
    
    % display
    disp(' ');
    disp(['map.totalNumberOfSites = ' num2str(totalNumberOfSites) ';']);
    disp(['map.totalNumberOfSpikes = ' num2str(totalNumberOfSpikes) ';']);
    disp(['map.spikesPerSite = ' num2str(spikesPerSite) ';']);
    disp(['map.normTotalNumberOfSpikes = ' num2str(normNumber) ';']);
    disp(['map.meanWeightedDistanceFromSoma = ' num2str(meanDist) ';']);

    handles.data.map.totalNumberOfSites = totalNumberOfSites;
    handles.data.map.totalNumberOfSpikes = totalNumberOfSpikes;
    handles.data.map.spikesPerSite = spikesPerSite;
    handles.data.map.normTotalNumberOfSpikes = normNumber;
    handles.data.map.meanWeightedDistanceFromSoma = meanDist;
end

% set up figure
if figFlag
    hfig = figure(109);
    set(hfig, 'Position', [65   695   560   420], 'DoubleBuffer', 'on', 'Color', 'w');
    pltrows = 2;
    pltcols = 2;
    pltnum = 0;
    colormap jet2;
    
    % plot the average map
    pltnum  = pltnum + 1;
    hplt = subplot(pltrows, pltcols, pltnum);
    himg = imagesc(mapAvg);
    set(himg, 'XData', [min(colVals) max(colVals)], 'YData', [min(rowVals) max(rowVals)]);
    addPixInterpOption(himg);
    xlim([min(colVals) max(colVals)]);
    ylim([min(rowVals) max(rowVals)]);
    pbaspect([1 1 1]);
    colorbar;
    hold on;
    hsoma = plot([0 0], [0 0], 'wo');
    ylabel('Distance (um)');
    title('R-click for options');
    
    % plot distances vs responses
    pltnum  = pltnum + 1;
    hplt = subplot(pltrows, pltcols, pltnum);
    cla;
    plot(distances, responses1, 'bo');
    hold on;
    
    herr = errorbar(Xuniq, Aavg, Asem, 'k-');
    try % modified for version compatibility GS20060530
        set(herr(2), 'Marker', 's', 'MarkerFaceColor', 'none', 'LineWidth', 2);
    catch
        set(herr, 'Marker', 's', 'MarkerFaceColor', 'none', 'LineWidth', 2);
    end
    xlabel('Distance (um)');
    ylabel('Number of spikes');
    title('Single cell resolution');
    
    % plot folded data, gaussian fit
    pltnum  = pltnum + 1;
    hplt = subplot(pltrows, pltcols, pltnum);
    
    [r1,c1] = size(responses2);
    responses2 = reshape(responses2, r1*c1, []);
    [r2,c2] = size(distances2);
    distances2 = repmat(distances2, c1/c2, 1);

    try
        out = responseDistanceGaussFit(hplt, distances2, responses2);
        fitting_parameters = out % for command-window display

        % how can you get the results in a usable data format?
        % i.e., how do you unpack the cfit class?
        % answer:
        coeffs = coeffValues(out);

        hfit = findobj(get(gca, 'Children'), 'Color', 'r');
        fitVals = get(hfit, 'YData');
        fitX = get(hfit, 'XData');
        peakVal = max(fitVals);
        minVal = min(fitVals);
        halfMax = (peakVal - minVal)/2;
        I = find(fitVals>halfMax);
        %     figure, plot(fitVals(I), 'ro')
        fullWidth = round( fitX(I(end)) - fitX(I(1)) );
        disp(['full width at half-max: ' num2str(fullWidth)]);
        disp(['half-width at half-max: ' num2str(fullWidth/2)]);

        handles.data.map.fwhm = fullWidth;
        handles.data.map.hwhm = fullWidth/2;
    catch
        disp(lasterr);
        handles.data.map.fwhm = [];
        handles.data.map.hwhm = [];
    end
end


%     % plot distances vs responses
%     pltnum  = pltnum + 1;
%     hplt = subplot(pltrows, pltcols, pltnum);
%     cla;
%         
%     % excitation resolution analysis (see also excitationResolution)
%     hax = gca;
%     rMax = max(get(gca, 'Xlim'));
%     r = 0:rMax;
%     
%     D = 50000; % neuronal density; units: neurons mm^-3
%     D = D*1000^-3; % convert mm^3 to um^3:
%     
%     step = 1;
%     r0 = 0:step:(rMax-1); % calculate the number of neurons as a function of r
%     r1 = r0 + 1*step;
%     dA = pi * (r1.^2 - r0.^2)/step;
%     z = 100;
%     N = dA * z * D;
%     k = max(N)/max(r1); % slope
%     
%     R = responses1.*(k*repmat(distances, 1, P));
%     plot(distances, R, 'bo');
%     hold on;
% 
%     hplot = plot(r1, N, 'r-');
%     
%     Y = k * U;
%     Ravg = Aavg.* Y;
%     Rsem = Asem.* Y;
%     herr = errorbar(Xuniq, Ravg, Rsem, 'k-');
%     set(herr(2), 'Marker', 's', 'MarkerFaceColor', 'none', 'LineWidth', 2);
%     xlabel('Distance (um)');
%     ylabel('Number of spikes');
%     title('Single flash resolution');
% %     Xval = find(R == max(R));
% %     disp(['X val of peak: ' num2str(max(Xval))]);
%     
% end

% X = distances;
% Y = responses1;
% % X = distances2;
% % Y = responses2;
% % --------------------------------------------------------------
% function out = getData(group)
% switch group
%         
%     case 'EP_adult_rat_L4'
%         list = [
%             'g0656P1'; ...
%                 'g0665P1'; ...
%                 'g0678P1'; ...
%                 'g0655P1'; ...
%                 'g0669P1'; ...
%                 'g0676P1'; ...
%             ];
%         data = loadDataFromList(list, 'apNumber');
%         out.maps = data.field;
%         vals = data.field
%         out.spacingXY = [100 100];
%         out.somaXYZ = [0 0 0];
%         
%     case 'EP_adult_rat_L5A'
%         list = [
%             'g0645P1'; ...
%                 'g0647P1'; ...
%                 'g0649P1'; ...
%                 'g0651P1'; ...
%                 'g0652P1'; ...
%                 'g0646P1'; ...
%                 'g0648P1'; ...
%                 'g0650P1'; ...
%                 'g0653P1'; ...
%                 'g0654P1'; ...
%         ];
%         data = loadDataFromList(list, 'apNumber');
%         out.maps = data.field;
%         vals = data.field
%         out.spacingXY = [100 100];
%         out.somaXYZ = [0 0 0];
%         
%     case 'EP_adult_rat_L4_092404'
%         list = [
%             'g0524P1_15mW'; ...
%                 'g0525P1_15mW'; ...
%                 'g0526P1_15mW'; ...
%                 'g0527P1_15mW'; ...
%                 'g0528P1_15mW'; ...
%             ];
%         data = loadDataFromList(list, 'apNumber');
%         out.maps = data.field;
%         vals = data.field
%         out.spacingXY = [25 25];
%         out.somaXYZ = [0 0 0];
%         
%         
% end
% 
